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Abstract 



In [3], we proposed a numerical regularized moment method of arbitrary order 
(abbreviated as NRxx method) for Boltzmann-BGK equation, which makes numerical 
^\ , simulation using very large number of moments possible. In this paper, we are further 

CNJ ' exploring the efficiency of NRrrx method with techniques including the 2nd order 

HLL flux with linear reconstruction to improve spatial accuracy, the RKC schemes 
' to relieve the time step length constraint by the regularization terms, and the revised 

O |. Strang splitting to calculate convective and diffusive terms only once without loss 

of accuracy. It is validated by the numerical results that the overall efficiency is 
significantly improved and the convergence order is kept well. 



1 Introduction 



J> ■ In 1949, the moment method was introduced by Grad [5] as a technique to approximate 

0^ , the Boltzmann equation in a macroscopic view. With proficient mathematical skills, he 

derived the famous 13-moment system, but the system is problematic: the hyperbolicity 
I/"") ■ can only be obtained in the neighbourhood of Maxwellian, and the structure of shock wave 

is non-smooth when the Mach number is large (see e.g. [13, 6]). During a long time, the 
moment method suffers lots of criticism, and very few progresses are made before 1990s. 
In the recent twenty years, a number of new ideas based on the moment method have 
come out, such as the Jin-Slemrod regularization of the Burnett equations [9], the COET 
(Consistently Ordered Extended Thermodynamics) method [12], the order of magnitude 
approach [15], and et al. In [17], regularization based on 1st order Chapman-Enskog 
expansion was considered to remedy the defects of Grad's 13 moment equations. The 
regularized system obtained therein was referred as R13 equations. Soon, the R20 and 
R26 equations are studied [11, 7] and the aspiration to extend this method for system 
with more moments [16] was called as Hxx by the authors of [17]. In [3], we proposed a 
numerical regularized moment method of arbitrary order for Boltzmann-BGK equation. 
The method is abbreviated as "NRxx method" later on for convenience following the 
tradition. The NRxx method makes it possible to investigate large moment systems, by 
direct numerical approximation of the regularized moment systems without deriving the 
explicit forms of the moment equations. The regularization in [3] is similar to the original 
derivation of the R13 equations [17]. And later in [4], the idea of COET [12] is adopted 
to revise the regularization terms, which results in a parabolic system. 
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Though one can escape from deriving the complex moment system, and the develop- 
ing of the simulation program is greatly simplified by the NRxx method, the numerical 
efficiency of the method in [3] should be further explored. It has been verified that the 
computational time is linear in the number of moments. However, a 1st order HLL nu- 
merical flux, which is over diffusive, was used in [3] for the transportation part, such that 
we have to use quite fine spatial grids to reduce the numerical error to a moderate level. 
Noting that the regularized moment system is parabolic, the time step is quadratic in the 
grid size. It turns out that in the ID case, the total computational cost is cubic in the 
number of spatial grids. With a quite fine spatial grid to achieve enough accuracy, the 
numerical simulation therein is rather time consuming. 

In this paper, we focus on improving the efficiency of the NRxx method in [3] with 
some state-of-art techniques available to us. Precisely, three numerical techniques are 
employed: 

1. Linear spatial reconstruction. The piecewise linear spatial reconstruction is able to 
provide a 2nd order HLL flux and greatly reduces the numerical diffusion. Owing to 
the absence of analytical expressions of the moment equations, the reconstruction 
needs to be done carefully. A conservative reconstruction is proposed, and it is 
numerically verified to have achieved high resolution. Although the final scheme is 
still of first order, the numerical error is greatly reduced. Thus much less grids are 
used in the computation. 

2. RKC time integration scheme. The RKC time integration scheme [21] is adopted 
to enlarge the time step sizes. The RKC method is a Runge-Kutta type method 
originally designed for diffusive PDEs to provide large stability regions, while the 
regularized moment equations are convective-diffusive problems, where imaginary 
parts appear in the eigenvalues of the right hand side of the semi-discrete system. 
As is well known, the RKC schemes contain a parameter called as the damping 
factor, which allows a small imaginary perturbation of eigenvalues and is usually 
selected as a small positive value. In order to make the RKC scheme compatible 
with the convection terms, we use a large damping factor in our scheme to ensure 
stability. The final time step size is equivalent to the grid size, and the number of 
internal time steps is equivalent to the square root of the grid number. Thus the 
total number of time steps is essentially reduced, and no instability phenomenon 
apprears in our numerical test. 

3. Revised Strang splitting method. With enlarged time step, the numerical accuracy 
can be harmed by convection-collision splitting. Therefore, the Strang splitting is 
utilized to win higher order of accuracy in the time direction. Usually, the implemen- 
tation of the Strang splitting requires twice calculations of the collision term in one 
time step. For this problem, we combine the two steps of collision in the successive 
time steps, so only once calculation of the collision term is performed in a time step. 

Several numerical examples are carried out to show the efficiency of our algorithm. Our 
prediction of the computational cost together with the order of convergence is validated 
by numerical experiments. We also compare the current method with the scheme with- 
out linear reconstruction to demonstrate much higher resolution for the one with linear 
reconstruction. 

The rest of this paper is arranged as follows: in Section 2, we give a brief review of the 
Boltzmann-BGK equation and the NRxx method. In Section 3, the linear reconstruction 
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is added to the HLL numerical flux. In Section 4, the time step size is enlarged by the 
RKC method, and the accuracy is improved by the Strang splitting method. In Section 
5, numerical examples are carried out to make illustrations. Some concluding remarks are 
given in Section 6. 

2 The NRn method 

2.1 The Boltzmann-BGK model 

In the kinetic theory, it is generally accepted that the Boltzmann equation is able 
to describe the fluids accurately. Due to the complexity of the collision term, in the 
computational field, several simplified collision operators are adopted, among which the 
BGK model [1] is the simplest but useful. The Boltzmann-BGK equation reads 

^+S-V x f = v(f-f M ), x,£eR D , (2.1) 

where / is the distribution function, v is the collision frequency, and /m is the local 
Maxwellian defined by 

, / ^ p(t,x) ( \£ — u(t, x)\ 2 \ 

Here p, u and are local macroscopic variables which represent the density, velocity and 
temperature respectively. They can be calculated from the distribution function / by 

p(t,x)= [ f(t,x,£)d£, (2.3) 

p(t,x)u(t,x) = f £f(t,x,£)d£, (2.4) 

p(t,x)u(t,x) + Dp(t,x)9(t,x)= ( \£\ 2 f(t,x,£)d£. (2.5) 

2.2 A first-order scheme for the NRxx method 

The NRxx method is raised in [3, 4] as a new tool for the computation of large moment 
systems. It is based on an Hermite expansion of the distribution function: 

f(t,x,£) = fa(t,x)n Bit>x) , a M-J^M , (2.6) 

where a is a D-dimensional multi- index. Hg^s and / Q 's act as the basis functions and 
the corresponding coefficients respectively, and 1-Le,a is defined by 

/ I |2\ D 

Hg, a (v) = (2tt)-t e~2 exp ( ) ]J He ad (v d ), (2.7) 

^ ' d=i 

where He n (x)'s are the orthogonal Hermite polynomials defined by 

He n {x) = (-l) n ^ exp(-x 2 /2). (2.8) 
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This idea originates from [5] where the 13-moment equations are derived. In order to get 
a finite system, the NRxx method [4] chooses an positive integer M and approximate f a 
with \a\ = M + 1 by 

Now that all / a 's with \a\ = M + 1 have been related with /^'s with |/3| = M, we can get 
a closed moment system of all / 7 's with | — y | ^ M by putting the expansion (2.6) into the 
BGK equation (2.1). 

The explicit expressions of such moment systems can be written in a uniform style for 
any choice of M, but those expressions are not convenient for computation, since they are 
not in a conservative form. In order to avoid an intricate process to obtain balance laws, 
the NRxx method treats the distribution function (2.6) as a whole, instead of considering 
each moment f a as an individual variable. For the construction of an applicable scheme, 
it is necessary to provide a method which is able to apply addition or subtraction on the 
following two distributions: 



= E f^ n ^ (^lP) ' = E KaUo^ (^Jp) 



(2.10) 



Additionally, due to the cutoff, only /i, a 's and /2,0's with \a\ ^ M + 1 are known. In [3], 
the authors proposed a homotopic method to calculate a new representation of f\: 

Hi) = E koMe 2 , a (^jp) , (2-11) 

where the coefficients /i jQ 's with |a| ^ M + 1 can be worked out by solving the ODE 
system 



(2.12) 



, D 

—F a = [1 - tR(t)] 2 E [R(r)0iF a . 2ed + (ui, d - u 2 ,d)VhMF a _ ed 
I ^(0) = /i, Q 

for all |a| ^ M + 1 until r = 1, and setting /i jQ , = F a (l). In Eq. (2.12), R(r) is defined 



as 



i?(r) = , V ^ - (2.13) 

The system (2.12) will be solved by a Runge-Kutta type scheme. Once the form (2.11) is 
obtained, f\ + f 2 or f\ — f 2 can be calculated naturally. 

By solving (2.12), we are able to represent / for any u' and 9' as 

/«)= E ^>(^r) + ---> ( 2 - 14 ) 

once there is one such representation known for some particular u' and 6'. Here the ellipsis 
means the remaining coefficients are unknown. For any u' and 0' and the associated 
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representation of / (2.14), we have 

P= [ /(£)d£ = /o 

PUd = 



o- 



/ e«*/(0^ = /o«d+4> d=i,---,D, 

D 

p\u\ 2 + Dp6 = |£| 2 /(£) d£ = p\uf + Y^ffd + + 2/£ 



(2.15) 



When u' = u and 0' = 0, (2.14) is called as the standard representation of /. Note that 
the coefficients in (2.6) and (2.9) are in the sense of standard representation. 

Another important technique is to calculate the flux £/ based on the Hermite expansion 
of /. Since only the moments with \a\ ^ M + 1 are known, £/ can only be accurately 
given upto the M-th order moments. Suppose / is presented as (2.14) for some u' and 6', 
and we let 



\ol\<M 



U 



+ ••• 



,D. 



(2.16) 



Then, by making use of the recurrence relation of the Hermite polynomials, we have 



F' 

3, a 



O'f' + u'j' a + (ay + 1)/' \a\ < M. 



(2.17) 



Here f ' . is taken to be zero if a = 0. 

Once the linear operations between discrete distributions and the calculation of the 
flux are applicable, the convection progress can be simulated by a Riemann-solver-free 
finite volume method, such as Lax-Friedrichs scheme and HLL scheme (see e.g. [10]). 
Suppose the problem is one-dimensional, and the spatial mesh is uniform with mesh size 
Ax. We denote by /" the distribution function on the i-th grid at the time step t n . Other 
symbols such as u™ and Of are defined similarly. Then, a first-order scheme with HLL 
numerical flux can be described as follows: 

1. Let n = and set /"(£) to be the initial value, which is in its standard representation. 

2. Apply (2.9) to obtain the (M + l)-st order moments: 



ft 



an tn _ tn 

°i Ji+l,a— ei J i— l,a— ei 

2Alc 



, \a\=M + l. 



(2.18) 



3. Solve the convection part with the HLL scheme: 



(2.19) 



where 



+1/2' 



f urn, oof H 
Af +1/2 ei/r(o-Af +1/2 6/f +1 (^ 

A i+l/2 A i+l/2 

Af +1/2 A^ 1/2 [/f +1 (g)-/f($)] ^ n< . X ij 



A i+l/2 A i+l/2 



0> A* 



+1/2- 



(2.20) 
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In Eq. (2.20), the signal speed A^ 1 / 2 and A^jwg are approximated by 



1/2 clliu 

\+l/2 - Ului \ 



Af, , /0 = min (n?, - C M+ iV0f, - C M +i,/^) , 

V / J ^ (2.21) 

A? 



2 - max {u^ + CW+iVflf, + C M+ i^f&^i) , 



where Cm+i is the maximal root of Hermite polynomial HeM+\{ x )- Eq- (2-21) is 
also used to determine the time step At. We refer to [3] for details. Note that (2.19) 
is only accurate up to the Mth order moments. 

4. Solve the collision part analytically: 

(a) For each /"*(£), get its standard representation using (2.12) and (2.15). We 
suppose the result is 

/no = E /Pv+sa \^-j=- l 

(b) Multiply all coefficients /"*, \a\ ^ 2 by exp(— uAt). The result is 

n+l 

/r +1 (o = E I *-^4t - i 

|a|<M 

where 

= I ^ H < 2 ' (2.24) 

5. Increase n by 1 and return to step 2. 

In step 3, since diffusion terms exist in the equations, the time step satisfies At = 0{Ax 2 ). 
The whole scheme is of first order. 



3 A high resolution scheme 



As is well known, the first-order HLL flux (2.20) adds excessive diffusion to the nu- 
merical solution in a general case. In order to reduce numerical diffusion, the technique 
of linear reconstruction is introduced below to the finite volume method. Suppose the 
boundary between the i-th and {i + l)-st cells is located at x = Xj+i/2- Our aim is to 
construct two distributions //j_ 1//2 (^) an d f^_i/ 2 (£) based on all cell averages fi to approxi- 
mate the left and right limit of f(x, £) at point x = x i+1 / 2 - In this section, the superscript 
n will be omitted since the reconstruction is only applied on the same time step. Thus 
the numerical flux in (2.19) can be re-formulated as 



Gi+l/2(£) 



-1/2 



(0, 



< Af , 



1/2' 



A i+ i/ 2 £i/i + i/2(£) ~ \+i/2£i/i+i/2(£) 



+ 



£i/£i /2 (0, 



-1/2^ 
A i+l/2 
^f+l/2^i+l/2L 



A i+l/2 
*R 



[^1/2(0-/^1/2(0] 



(3.1) 



A i+l/2 



A i+l/2 



X i+l/2 < < A i+l/2' 



^ A* 
j+i 



/2' 
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where 

\+l/2 

^+1/2 - luaA V"l,i+l/2 



A f+l/2 = min { n M+i/2 - C M+lyj0f +1/2 , ""£+1/2 - ^+1^+1/2} ' g 2 ^ 

A^ii /o = max (uf „■_,_-, /0 + Cm+i n M+i/2 + Cm+i ^^-1/2} ■ 



During the reconstruction, different methods are applied to the convection part (\a\ ^ 
M) and the diffusion part (\a\ = M + 1). 

3.1 Reconstruction for the convection part 

For the convection part, it is important that the quantities used in reconstruction 
are conservative variables. Before reconstruction, according to the algorithm described in 
section 2.2, we already have the standard representations for all distributions /«(£), and 
the coefficients are assumed to be fi )OC , \a\ ^ M. The simplest idea is to use fi )OC together 
with ui and Q{ to make linear reconstruction: 

fi—l/2,a = fi,a ~ 9i,aAx/2, fi+\/2,a = 9i,aAx/2, 

u f-i/2 = u i ~ 9 Ax/2, uf +1/2 = Ui+ g i Ax/2, (3.3) 



£1/2 = °i- df +1/2 = 9 t + g t Ax/2, 



and 



f^-1/2^) — E ^l/2,a% 1/2 ,a 




+ ••• 



fi+1/2^) ~ E /il/2,a«^ 1/al a ( fpZ~~^ ) " 

\ A/ i+l/2 



(3.4) 



where and ^ are constants, and gfj is a constant vector. However, this method leads 
to incorrect numerical results since none of the variables in Eq. (3.3) is conservative. 

We take the original idea of the NRxx method, and consider the distribution function 
as a whole. Thus, linear reconstruction means taking the following approximation of 
/* 1/2 (€) and /£ 1/2 (£): 

/£i/2«) = " /ii/a«) = /*«) + (3-5) 

Here <?j(£) is a distribution. Obviously, the most convenient representation of <?i(£) is 

= E ( J + • • • ■ (3-6) 



Thus no ODEs are to be solved during the calculation of (3.5). Now, the coefficients ff^a's 
can be naturally given if we represent and as 



|a|<M 



(3.7) 



|a|<M 
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In the implementation, we use the simplest minmod slope limiter for reconstruction 

, I fi+l,a ~ fi,a fi,& ~ fi-l,a I /„ \ 

g ha = minmod j — , ^ j . (3.8) 

This reconstruction is a conservative reconstruction, since 

[ Xl+1/2 fi(x,$) = Axfi(Z), VzeZ, (3.9) 

where fi(x,£) is a linear function of x defined on (^-1/2)^1+1/2) as 

fiM) = fi(Z) + &(0 (z - Xt ' 1/2 ^ Xm/2 ) , x e ( Xi _ 1/2 ,x i+1/2 ). (3.10) 
However, this condition is not satisfied by the reconstruction (3.3). 

3.2 Reconstruction for the diffusion part 

The diffusion terms (2.9) provide approximation to the (M + l)-st order moments. 
Since (2.9) is in the sense of standard representation, we first need to get the standard 
representations of //Li/ 2 (0 an< ^ /i+1/2 an d we use (3.4) to denote the results. The 
reconstruction of 7/^/2 a an d /i+1/2 a with \a\ = M + 1, which are involved in the ellipsis 
of Eq. (3.4), is a direct discretization of (2.9): 

rR ^i— 1/2 fi,ct—ei fi—l,a—ei 

e l L 1/2 _ H=M + 1. (3.11) 

j.^ i+l/2 /j+l,a— ei Ji,a—ei 

h+l/2,a ~ ~~L Xx ' 

^+1/2 ^ X 

Now we give a general discussion on the order of accuracy. The time splitting introduces 
an error of mag nitude O(At) = 0(Ax 2 ). With (3.5), (3.6) and (3.8), the numerical flux 
(3.1) turns out to be a second order numerical flux. However, as discussed in [18, 20], due 
to the one-sided approximation of the diffusive gradients (3.11), the final accuracy only 
appears to be the first order. However, comparing with the original scheme described in 
Section 2.2, numerical error is significantly reduced by the linear reconstruction. 

4 Enlarging the time step 

The regularization of the moment method introduces diffusion terms into the system, 
which yields a relatively small time step At = 0(Ax 2 ). In order to enlarge the time step 
length, we use the RKC time stepping in the temporal discretization. In this section, a 
large time-stepping scheme with 2nd-order time integration will be proposed. 

4.1 The RKC time-stepping 

The RKC method is a series of explicit Runge-Kutta schemes for parabolic problems 
with large stability region and good internal stability. Our aim is to use At = O(Ax) in 
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our algorithm without producing errors larger than the original method. Thus a second- 
order RKC scheme is needed. The s-stage second-order RKC formula for the ODE system 



w'(t) = F(w(t)) (4.1) 

was deduced in [21] as 
Wo = w n , 

W l = W + /iiAtFo, 

Wj = (1 - Hj - Vj)W + VjWj-x + + fij&tFj-! + 7jAtF , j = 2, ■ ■ ■ , s, 

w n+l = W s , 

where Fk = F{Wk), fix = b\LJi, and 

Mj = *r , ^' = -7 , Mi = -r , Ti = -o-j-iHj, J = 2, • • • , s. (4.3) 

: 6 i-2 

Here the parameters dj, bj, ljo, oj\ are relevant to a manually selected damping factor e. 
They are given by 

uo = l + es 2 , u 1 = r a (ui )/T?(uio), 

b = h = b 2 , bj = TjVo)/(T> )) 2 , j = 2, • • • , s, (4.4) 
(ij = l-bjTj(uo), j = 1, ••• ,s- 1, 

where Tj{x) is the first kind Chebyshev polynomials 

Tj(x) = cos(j arccosx) = cosh(j arccoshx), j = 0, 1, 2, • • • . (4-5) 

which can also be defined by the recurrence relation 

T Or) = 1, T 1 (x) = x, T j+1 (x) = 2xT j (x) - T^i(x), j = 1, 2, • • • . (4.6) 

The damping factor e is often chosen as a small positive value such that a small imaginary 
perturbation of the eigenvalues of F'(w) is allowable. In [8], the authors suggest that e be 
chosen as 2/13, which results in a reduction in the stability boundary of about 2%. 

For advection-diffusion problems, the imaginary part of the eigenvalues of F'(w) may 
be large. Some analysis of the eigenvalue structure of upwind finite volume methods can 
be found in [19], where the authors show that the imaginary parts of the eigenvalues are 
less than 2.0 for the second-order scheme if the time step satisfies the CFL condition of 
a pure advection problem. In order to ensure the stability, we follow [22] and use a large 
damping factor e = 10. Thus the stability boundary is approximately given by 

/3(s) w 0.34(s 2 - 1). (4.7) 

In our implementation, the RKC method is only applied to the finite volume scheme 
(2.19). The time step is determined by 

< CFL, (4.8) 

Ax 

where A max is defined as 

A max = max {Af +1/2 , Ag_ 1/2 } . (4.9) 
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Then, we find the smallest positive integer s satisfying 

and use the s-stage RKC method instead of the forward Euler scheme in (2.19). The 
inequalities (4.8) and (4.10) give At = O(Ax) and s = O (l/^/Axj , which leads to the 
following estimation of the total computational time: 

t com ^-!—-^ = 0(Ax- 2 - 5 ), (4.11) 
At/s Ax 

where L is the length of the computational domain, and T is the finishing time of the 
computation. 

4.2 The Strang splitting 

When the RKC time stepping is used, according to (4.8), the time step has the same 
magnitude as the grid size. Thus, when the gas is dense, we can still have only the 
first-order accuracy due to the convection-collision splitting, which introduces an error 
O(At) = O(Ax). In order to restore the time integration to a second-order one, the 
Strang splitting technique [14] is employed. In our algorithm, a direct usage of the Strang 
splitting can be described as follows: 

1. Let n = 0. 

2. Determine the time step At n . 

3. Solve the collision part over a half time step of length At n /2. 

4. Solve the convection part using the second-order RKC scheme over a time step of 
length At n . 

5. Solve the collision part again over a half time step of length At n /2. 

6. Increase n by 1 and return to Step 2. 

This scheme requires twice calculation of the collision operator in one time step, but 
decrease the time integration error to 0(At 2 ). Here we introduce a method equivalent to 
the above Strang splitting scheme, but only once calculation of collision term is needed 
per time step. 

According to (4.8) and (2.21), the time step determined in Step 2 is only relevant to u 
and 9. Since the collision steps (Step 3 and Step 5) do not change these two quantities, we 
can actually calculate At n+ \ before Step 5. Thus the above algorithm can be rearranged 
as follows: 

1. Let n = and determine the time step At®. 

2. Solve the collision part over a half time step of length At n /2. 

3. Solve the convection part using the second-order RKC scheme over a time step of 
length At n . 

4. Determine the next time step At n +\. 

5. Solve the collision part over a half time step of length At n /2, and then solve the 
collision part again over a half time step of length Ai n +i/2. 
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6. Increase n by 1 and return to Step 3. 

Recall that for the BGK model, the collision part can be solved analytically as (2.23) and 
(2.24). Therefore, we can replace Step 5 by 

5. Solve the collision part over a time step of length (At n + Ai n+ i)/2. 

Thus, the numerical result is identical to the original Strang splitting scheme, but only 
one collision step is performed in a complete time step. 

5 Numerical examples 

In this section, two numerical examples are presented to validate the efficiency and 
accuracy of our algorithm. In all the tests, the collision frequency v is given by a simple 
form p/ ' Kn, which is corresponding to the Maxwell molecules. Here Kn is the global 
Knudsen number, and it is slightly different from the Knudsen number defined in [2] 
(denoted by Kn') by 

Kn = -a -Kn. (5.1) 

5 V 7T 

The CFL number is chosen to be 0.95. All the computations are performed on the Dell 
OptiPlex 755 desktop computer with a dual-core processor and CPU speed 2.33GHz. 

5.1 An example with smooth solution 

The first example is a repetition of the ID periodic problem in [18, 3]. The computa- 
tional domain is [—1, 1] and the boundary condition is assumed to be periodic. The initial 
condition is 

p (x) = 2 + i cos(7rx), u Q (x) = ^1 + ^sin(7nc), isin(7rx),0^ , p (x) = 1, (5.2) 

and the distribution is the local Maxwellian everywhere. The computation is stopped at 
t = 0.4. We set Kn = 0.5 and M = 3,6, 9. The numerical results for the density p and 
temperature 9 are plotted in Figure 1. Since the Knudsen number is large, the profiles of 
density and temperature differ a lot for different moment systems. 

In order to examine the efficiency of our algorithm, we discretize the problem on a 
series of spatial grids with grid numbers ranging from 10 to 200. Using N = 2/ Ax, the 
analysis at the end of Section 4.1 gives 

t com = 0(N 2 - 5 ), At = 0(l/N), At/s = 0(l/N 1 - 5 ). (5.3) 

These results are validated by Figure 2. We can see that the large time step sizes are 
achieved, and the results are still stable. This leads to a remarkable reduction of the 
computational time. 

Also, the first-order convergence rate is illustrated in Figure 3, where the "exact" 
solution is obtained on a mesh with 800 grids and the L 1 errors are shown. We would like 
to remark that although the scheme is still of first order, the magnitude of error is much 
smaller than the original method. In [3], where the HLL scheme without reconstruction 
was utilized, the authors used 1000 grids to obtain numerical results with similar resolution 
as those in Figure 1. Note that small grid size implies smaller time steps, which leads to 
a huge computational cost. 
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(e) M = 9, density (f) M = 9, temperature 

Figure 1: Numerical solutions of Problem 5.1 using 200 spatial grids 
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Figure 2: Computational costs and time step sizes for Problem 5.1 on different spatial 
grids. The x-axis is the logarithm of the grid number N. 
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Figure 3: Convergence rates for different moment systems 
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5.2 The shock-tube test 



In this example, we show that our method is able to achieve high resolution when sharp 
layers exist in the numerical solution. Here a Riemann shock-tube problem is considered. 
It has been studied by Yang and Huang in [23] using the discrete ordinate method. The 
initial states are 

( P ,u,e) = \ fa"*'*)' x<0 n l> (5.4) 

^ ' \ (p r ,u r ,e r ), X > 0.5 v ; 

with _ 

pi = 0.445, u t = (0.698\/2, 0, 0) T , t = 13.21, 

rp (5.5) 

Pr = 0.5, u r = (0, 0, 0) T , 6 r = 1.9. 

The Knudsen number are selected by setting Kn' = 0.001 in (5.1). The computational 
domain is set to be [0,1], and we solve the problem until t = 0.1314/\/2 ~ 0.09291 in 
order to match the results in [23]. 

Since Kn is small, only the case M = 3 is considered here. If large M is used, the 
results are nearly identical to the current case. Some results are listed in the left column 
of Figure 4, whose validity can be confirmed by comparing them with those in [23]. In 
order to see the effects of reconstruction, we set g^ a = in (3.8) and rerun the program. 
The results are in the right column of Figure 4. It is obvious that the left column provides 
much higher resolution near the shock wave, while the right column is even unable to 
achieve the correct peak value for TV = 100 and N = 200. 

6 Concluding remarks 

An efficient numerical scheme with high resolution for the NRxx method has been 
presented. Since the NRxx method gives a convective-diffusive system, we not only per- 
form the linear reconstruction to gain a high spatial resolution, but also use the RKC 
schemes and the Strang splitting method to enlarge the time step while maintaining the 
order of accuracy in Ax. In the future work, we are extending it into the 2D case with 
unstructured grids, together with the specularly reflective boundary conditions. 
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